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Abstract 


We develop and test a 1— point closure turbulence model with the following features. 

1) we include the salinity field and derive the expression for the vertical turbulent 

diffusivities of momentum K . heat K, and salt K as a function of two stability 

m’ h s 

parameters: the Richardson number Ri (stratification vs. shear) and the Turner number R^ 
(salinity gradient vs. temperature gradient). 

2) to describe turbulent mixing below the mixed layer (ML), all previous models have 
adopted three adjustable "background diffusivities" for momentum, heat and salt We 
propose a model that avoids such adjustable diffusivities We assume that below the ML, the 
three diffusivities have the same functional dependence on Ri and R as derived from the 
turbulence model However, in order to compute Ri below the ML, we use data of vertical 
shear due to wave-breaking measured by Gargett et al. (1981). The procedure frees the 
model from adjustable background diffusivities and indeed we employ the same model 
throughout the entire vertical extent of the ocean. 

3) in the local model, the turbulent diffusivities K , are given as analytical functions of 

Ri and R„ 

P 

5) the model is used in an O-GCM and several results are presented to exhibit the effect of 
double diffusion processes. 

6) the code is available upon request. 
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I. Introduction 

For sake of completeness, we recall that the 0— GCM solve the dynamic equations for 
the mean velocity th, mean temperature T and mean salinity S: 

a u i + m. ( u i u j + 5Fj ) = - < la > 

+ < lb > 
f + ^(U iS + ^) = ... (ic) 

The velocity, temperature and salinity fields have also fluctuating components u-', T" and 
s" which produce the correlations uFut 1 (Reynolds stresses), uVT"(heat fluxes) and 
u-'s" (salinity fluxes). The challenge then is to construct such correlations so as to solve 
Eqs (la-c). To fix the ideas, we further write: 

(l d) 

(le) 

i 

dS 






(If) 


where S, =i(U ; -+U •) is the mean shear. The K , „ are the turbulent diffusivities for 
ij LJ hi m,h,s 

momentum, heat and salt. As discussed in paper I, they have the general functional form. 

K 2 C 


^m,h,s 2 f ^m,h,s 


(lg) 


where K and e are the turbulent kinetic energy and its rate of dissipation which in principle 

are given by two dynamic equations (the K-e model). The dimensionless structure 

functions S , must differ from one another so that: 
rxi }ii ju 


K 

m h r s 


(lh) 


In general we can write. 

S m,h,s = S m,h,s ( ' U ’ “s VS > <“> 

where a and o are the volume expansion coefficients a=-p A dp/ffT and a =p A dp/dS and 

I b lb 

where the shear VU can be generated either by external sources like in the mixed layer ML 
or by internal wave-breaking processes below the ML. If one introduces the Turner number 
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(2a) 


and the Richardson number Ri: 


VS“s f <*■*&> Ri = N h/ N S 


flT\. 


where 


N 


d T 


J h = s«TW- N 5 = 2 (W 
N2 = -?l = S“rf-8°sI= N h( 1 -V < 2b > 

we can rewrite (li) more concisely as 

S m,h,s= S m,h, S <V Ri > < 2c > 

(Clearly, we could have also defined Ri in terms of N 2 rather than just the thermal 
gradient. We have chosen the latter for reasons of presentation of the results). One must 
distinguish the following four cases: 

SF ( salt fingers, salty-warm over fresh-cold ): 


dS >Q 

3£ >0 ’ 


0T >n 

W >0 ’ 


R^>0, Ri>0 

R^<1 Stable, N 2 >0, R^>1 Unstable, N 2 <0 


(2d) 


DC ( diffusive convection, fresh-cold over salty-warm): 

Q 

3F <0 > 


^ <0) 


R^ >0, Ri<0 

R^<1 Unstable, N 2 <0, R^>1 Stable, N 2 >0 


(2e) 


(2f) 


DS ( doubly stable , fresh-warm over salty-cold) 

«* 0 5 T >0 

^• <0 > W >0 ' 

R^<0, Ri>0, N 2 >0, Stable 
DU ( doubly unstable, salty-cold over fresh- warm): 

dS n *r <a 
3i >0 ’ ^-<0, 

R^<0, Ri<0, N 2 <0, Unstable 
The stability/instability is predicated on the Brunt-Vaisala frequency N with N 2 >0 
(stable) and N 2 <0 (unstable). 

The general problem is to construct (2c) so as to encompass all four cases (2d-g) 


(2g) 
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First, there is ample evidence from laboratory and oceanic field data that show that is 
different from K g . In the SF case, the ratio K g /K^>l (Hamilton et al., 1989) while in the 
DC case, Kj i /K g >l (Kelley, 1984). Schmitt (1981) has also shown that the observed T-S 
relationship is not consistent with K^=K g . For a discussion and review of the importance 
of these processes and their extent in different parts of the ocean, see Turner (1967, 1973, 
1985), Schmitt (1994) and Zhang et al. (1998). In spite of this evidence, almost all O-GCM 
still assume 

K=K h (3a) 

Recently, attempts have been made to overcome (3a) but the task is not easy The mam 
difficulty is that in the absence of a model capable of encompassing all four cases, SF, DC, 
DS and DU, the only alternative is to employ laboratory and ocean data to build the 
functional form of the diffusivities to be used in an O-GCM. This is the approach 
employed by Large et al. (1994), Zhang et al. (1998, ZSH) and Merryfield et al. (1999, 
MHG) who used relations by Schmitt (1981) and Kelley (1984, 1990), among others. 

There is, however, an internal limitation to such a procedure since the available data 
refer to SF and DC but not to DS and DU which are also important (Duffy and Caldera, 
1999) Thus, away from the regions where SF and DC are active, the above authors take 

K m,h,s< DS - DU > =° < 3b > 

or, more precisely, they use a background diffusivity which is chosen primarily on grounds 

of numerical stability but whose physical origin must be an internal-wave breaking 

phenomenon This is clearly not a satisfactory situation especially in view of the fact that 

Since the studies by ZHS and MHG have shown the importance of double diffusion, the 

above procedure is certainly far better than (3a) but still not fully satisfactory. The goal of 

this paper is to consider the same problem but with a different methodology. 

We develop a turbulence model to compute the three diffusivities for momentum, heat 
and salt and construct the functions (lg) and (li) for the four processes SF, DC, DS and 
DU in the presence of an arbitrary shear. The inclusion of shear is quite relevant since is 
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known to hamper the SF mechanism (Linden, 1971, 1974a, b; Kunze, 1990) and yet, the 
above procedures do not account for shear since they expressed and K g in terms of only 
one stability parameter R^, rather than and the Richardson number Ri 

We present three models- 1) K and e are solutions of two dynamical equations (K-t 
model), 2) only one of them satisfies a differential equation while the other is taken to be 
the local limit of its dynamic equation and 3) both K and e are taken as the local limit of 
their respective dynamic equations. As we shall show, in model 3) all the relations are 
algebraic and one must solve a cubic equation. The numerical results correspond to 3). 

The structure of the paper is as follows. In II— VII we derive the general non-local, 
dynamic equations for the mean fields as well as the turbulent variables. In VIII we derive 
the analytic expressions for the turbulent diffusivities with only two non-local variables, 
the turbulent kinetic energy K and its rate of dissipation, e In IX-X we study the case of 
double diffusion without shear and show that the predictions of the model are in agreement 
with several laboratory and ocean data. In XI we give the complete analytic solution for 
the local model: we derive the algebraic expressions for the momentum, heat and salt 
diffusivities in the presence of arbitrary shear. In XII, we discuss the time scales. In XIII, 
we display several solutions of the model, specifically we plot the diffusivities K's and their 
ratios as a function of the two stability parameters, the Richardson number and the Turner 
number. In XIV-XVI we present the results of an 0— GCM with the above model where we 
use the same turbulence model below the mixed layer where the shear is no longer due to 
the external wind forcing but to a wave breaking mechanism. In XVII we present some 
concluding remarks. 

II. Continuity equation 

Following the formalism presented elsewhere (Canuto, 1997), the total velocity, 
density and pressure fields are split into mean and fluctuating parts as follows: 

u i = U i+ uV, p=p+p\ P=P+p\ p’=p'=0 (5a) 
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the Reynolds average <uj>=0. The relation between the two is discussed in 
(1997a). Using the equation for the density p 



we obtain, upon mass averaging, 


d - d d 
cTr 3t +u iac. 


Dt? + 4 yi=°’ 


D -d ii 5 

Dt = m + U i3x; 


III. Velocity Field. Mean and Turbulent Variables 


Consider the Navier— Stokes equations 

It " u i + 

J J 




j 


ij 


where a-j is the viscous stress tensor (v is the kinematic viscosity) 

_ ,d , d \ 2 <- d 

a ij - 

J 1 K 

Mass averaging (5e), we obtain the dynamic equation for the large scale flow U, 

+ r y) 


'’StV ~lx j (P,5 ii + r ij>-«5i 


where are the turbulent Reynolds stresses 


r- ■ = ou'.'u'-' = nR.. 
ij H i J H ij 

The kinetic energy of the large scale field 

K„ = iUjUj 

satisfies the equation (a,j= da/dx^; a^ 

^«-'-U 1 (P i i + r lu+ ^) 

The Reynolds stresses R.j satisfy the non-local dynamic equation (Canuto 1997): 

?( Et R ij + D ij> = A ij + B ij - J ij + *ij PD - ? e ij 
where the non-local term D- represents the flux of Reynolds stresses R-^: 

D u = r' i k l? R ijk + i*ijP' 5 k - Wj - Vi 1 

R ijk 

In Eq.(6a), the source term due to shear is represented by: 


Canuto 

(5c) 

(5d) 

(5e) 

(5f) 

(5g) 

(5h) 

(50 

(5j) 

(6a) 

(6b) 

(6c) 
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-V^k + W 

while the source (sink) term due to stratification is represented by 

;®ii = (?^ ik+ 7ip jk) p ik 

The fluctuating pressure p' gives rise to the pressure-velocity correlation 


7T:: = !!•: — o^-Ili 


n ij“ u i p ,J + u j p ,'’ ”ij = “ij “ J"ij" kk 
Finally, compressibility introduces a pressure-dilatation term 


(6d) 

(6e) 

( 61 ) 

(6g) 


where d=u'-' • is the "dilatation", while e-. is the dissipation tensor which we assume 
diagonal for its largest contribution originates in the small scales region: 


2 - 


(6h) 


£ ij - 

Below, we present the dynamic equation for e The trace of (6a) yields the equation for the 


turbulent kinetic energy K, 

K = ip* 1 pu'j'uv = iR n 

Ht K + D f = “ R ij U i,j + p' 2 P ,i + ~P A - f 
where D^K) is the non-local transport of K: 


D f E r' Lm? R kk> + - Vj 1 


(7a) 

(7b) 

(7c) 


IV. Concentration Equations 


Consider a model with two fluids of density pc and /?( 1 — c), where c is the 
concentration and p the total density of the fluid which satisfies (5c). The equation satisfied 
by pc is given by (no external sources) 



lt pc 

II 

O 

+ 

(8a) 

or alternatively, 






4t - ^ J i),i 

(8b) 

where is the diffusion flux 






T dc 

J i “ x cM 

(8c) 
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where * c is the molecular diffusivity of the c-field. Mass averaging, we derive 

pc = pC, pcUj = pCUj+ p^j (9a) 

where C is mean concentration and is turbulent concentration flux 

Cec, pul'c" = p<F (9b) 

Taking the mass average of Eq. (8b), we obtain the equation for the mean concentration C: 

< ]o *> 

To provide the turbulent flux 3^, we need a turbulence model. 


(ID 


V. Temperature field 

We begin with the equation for the entropy S (Landau and Lifshitz 1987), 

pT ft = + J il^ + a ij^ 

where /t is the chemical potential and 

q i = F i " J j(^ _T 1^1 P , C ) E F i " hJ i < 12 > 

r r* 

where Fj is the thermal flux. In the absence of diffusion, q^= F ( but here q } depends also on 
the gradients of the concentration as well as on the gradient of p. Using 


dS _ <9S , dT^dS, dc.dS, dp 


3T c7T'c,p 3T + dc 'T,p 3t + $p'c 


and the relations: 


rp ^S | dS | dp | 1 ( 8 1 dp | 

1 cTT ' c,p - p’ c?(VT,p - oT' ’ p Up'c,T ~ oT'p,( 

p,c 


(13) 

(14) 


one can transform Eq (11) as: 

yft^ T + fe> u i T ^ = ~ fxj F i + ^ + px c ^ k H k 

where the dimensionless function wis defined by: 

u = ~ E a T T 

Next, we take the mass average of (15a). Making use of the results derived in Canuto 
(1997) and recalling that a--u- • = pc, we obtain 

^pSr = - ( F i +F r "pXb + "Dt + “T p ,i - ^ + h + x c 


(15a) 


(15b) 
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(16a) 


where is the heat flux 

F i = y 7 ^ = P H j ( 16b ) 

The uF term in (16a) can be written using the second of (5b) and (22e) below. Eq (16a) is 

the generalized Bernoulli's equation with turbulence, diffusion and heat flux. When dealing 
with incompressible ocean turbulence, one can neglect the third-order term p'uV as being 
smaller than the second-order terms to which is summed; the last term can also be 
justifiably neglected as a kinematic term smaller than the remaining turbulence terms; the 
pTT term can also be neglected since d=<?Uj/<9xj=0 in an incompressible flow; the term pt 
cannot in principle be neglected since it represents the energy gained by the temperature 
field that is lost by the kinetic energy because of friction; thus, its presence is a 
consequence of energy conservation Thus, we have: 


pc plir ~ “ ( F ? +F i),i + ^Dt + ”iEE.) P + pt ( 16c ) 

Since D/Dt -dj d t +U| djdx^ one can probably neglect the uF term since it summed to the 
mean flow velocity lb which is expected to be larger. This further reduces (16c) to: 

P c pIJT = " lx^ F i +F i) + ^t + ' pt ( 16d ) 


In all O-GCM, this equation is however further reduced to: 



(16d) 


VI. Turbulence 

As already discussed, we need to evaluate the following second-order moments (2c). 
The equation for the first of them has already been given by Eq.(6a). As for the correlation 

<f> = fpc 112 = p$ (17) 

we first note that using Eq.(8a), pc 2 satisfies the dynamic equation (J- -EdpJ./dx^) 

’ <> 8a ) 
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Taking into account the relations 


p c 2 = p C 2 + pc 71 
/9UjC 2 = pUjC 2 + U^pc 172 4- pu^'c" 2 4- 2CpiT r c 7r 
the mass average of (18a) gives: 

It ^ + D f^) = ” ^ U i,i " ^i C ,i + ^,r CJ i,i 


Df(^) - (puj'c 02 ) 


Introducing the new function $ 


^pc 772 = p$ 


Eq (18d) simplifies to: 

?lDt $ + D f(*)l = -^i c , 1 + ^,r CJ i,i 

where the non-local transport of $ is given by 

D fW = 1 

The last two terms in (19b) will be evaluated as follows: 


(18b) 

(18c) 

(18d) 

(18e) 

(19a) 

(19b) 

(19c) 


UT - CX . = CX . + c"X . - CX . = cT ■ = c"XX . 

1,1 1,1 1,1 1,1 1,1 1,1 1,1 


(19d) 


which, using Eq (8c) and considering the incompressibility of the flow, becomes (for ease of 
notation, we employ <..> instead of an overbar) 

J( <9 2 c". 


«c <c ”3?j > 


(19c) 


Using a mathematical identity, we also have 

^c <c "lk? > = ~~ px z <( ^i) 2> + ^clkf ( 19f ) 

The last term represents the viscous-diffusion of (17), which we may consider small while 
the first term is the viscous-dissipation which we cannot neglect. We recall that in the case 
of momentum, the dissipation e, see Eq.(7b), 


< 19 e) 

is of the same form as the first term in (19f) and since one takes c=2K/r, by analogy we 
write 
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where r is correlation time scale to be discussed below Thus, finally, 

g t $ + D f ($) = -$ i C i -2r c 1 $ (20) 

Next, we consider the third term in (2c). Multiply (8a) by Uj, (5e) by c and add the 
results We obtain. 


ft^ cu i) + i^.(^ cu i u j) ~ F i c + u i J k,k 
J J 


where 


Recalling that 


d 


p i*-s.-«i + fe' i 


ij 


pcu- = p\]C + fmyr = pl^C + 


= j^UjC + Cr n + + 6^) + 

substitution into the mass averaged form of (21a) gives, after several steps, 

+ Vj-*j w u + ?1a i 

where the function A- is given by 

A i $r ?- F i C + “Pk,k 

After some algebra, we have: 

FT - FjC = - gpAjC 71- - <c"^ > + F j (vis)c - Fj(vis) C 


(21a) 

(21b) 

(21c) 

(21d) 

(21e) 

(21f) 

(22a) 


where F (vis) is the last term in (21b) and the dimensionless function Aj is given By 

A , = H p p -'f , H p = P (6?)-‘ < 22b > 

After some steps, we have: 


_ a2„M 

Fj(vis)c - Fj(vis) C = vp<c"^-^i > 


Since by definition 


p c IT = - p'c" 


(22c) 


(22d) 


use of the expansion 



gives 


( = - <T" + a c" 
P 


c" _ a x"c" - a^ 1 


(22e) 


(22f) 


The c 1 ' 2 term will be approximated with 2$ given by Eq (19a) while T"c" will be computed 
later. Thus, we have: 


A. = - gpA.(aT"c" - 2 a„ $) - <c M ^§ > 4- vp< c"— u i> + u-'J, ^ 

1 1 C <7X. 2 1 


5xj 


(22g) 


To evaluate the last term in (22g), we use (8c) and obtain 


"Fia = x c <(,u 


i, d 2 c ", 

'H' 


(24a) 


(24b) 


since by definition pu^=0 The last two terms in (22g) are therefore 

i^p<c"— — u i> + x3< u" ^“> 

which, using mathematical identities, we rewrite as 

*P(H-X c )$i jkk ~~P{^X Q )<%^ $%> + ^(^ c )[^ k (c"^0 -^ k ( u i^ k )] 

The first term represents the diffusion of We shall neglect the last term since one can 
argue that c" and the velocity gradient peak at different wavenumbers and there is 
therefore little overlap. As for the second term, it has a structure intermediate between 


(19g) and (19b) and it will therefore be written as 

?(^+X c )<^i J£ k > = ^ c r uc*i 

The last term we must compute is the pressure correlation term 

n? = <c"|£"> 

Using the analogy with the temperature case, we write 


n ?=^pi*i 


(24c) 


(24d) 


(25e) 


Finally, the complete equation for $• is 


D 


Ut*i+ D f (« , i)= - R.jCj- tjUjj - gA i (aTV T -2a c 4>) - + *x c 0+>'/X c )* ljkk 
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( 26 ) 


where we have absorbed r uc into r . Next, we consider the fourth function in (2c) which 
we generalize to 


= \^T T2 = pty 


(27a) 


First, we recall that, except for the last term, the temperature equation (16a) can be 
treated as in Canuto (1997), where the equation for rp is given by Eq.(26f). Thus, we must 
add the last term in (16a), 



(27b) 


in the derivation, one encounters the term 

<T "l^y" > = - < (w'b 2> + ^2 (27c) 

While the last term represents the diffusion of the potential energy ^T 772 , the first term 
represents the dissipation of it and we shall write it in analogy to the dissipation of 
turbulent kinetic energy, Eq.(19g), with a time scale to be discussed later. The final 
form of the dynamic equation for $ is: 


UT + p ' iD f ~ ~ c p lH i T ,i " 2r ^ + ,kk ~ T6*c Rc p^ Sx k 2& k 
where the heat flux Hj is defined in Eq.(16b). Next, we consider the second term in (2c), 
namely the heat flux (16b) Here too, the relevant dynamic equation was already derived in 
Canuto (1997), Eq (24a), to which we must add the last term of (16a) In the derivation 
one encounters a term analogous to (24b), specifically, 


25. 


dC <9T 2 


(27d) 


(27e) 


kT"— D i> + *<< *T> 

H 'H 

which we treat in a similar fashion. The final result is: 

D, H i+‘ P D f< H i> = - c P R ij T ,r H j u i.i - v A i T " - r P^i 

. ,/ , , \tt 25 p TT dC 5T 

+ “ ~8 R ^c U i ^3T k ^T k 

where the pressure term give rises to the relaxation term r^. Finally, we use the fact that 


(27f) 
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p Y T = -p r T T 


and the expansion (22e) to obtain 


T Tr = qT 72 - a c"T" 


(28a) 


(28b) 


so that (27f) becomes: 


Ct H i +<: p D f( H ,> = - c p R ij T r H j u i,r c p gAi( 2 a *-a c ?T") 
- r p0 H i + ■(‘ ,+ ^)H iikk - lS R Xc U i H k fl k 


(28c) 


Finally, let us consider the last term in (2c), the correlation between T" and c" We 


recall that in general 


and thus from Eq. (8b) 


dc DC , Dc" , n,d C , dc\ 

+u i^w i + w i ) 


(29a) 


(29c) 


(29d) 


^Sr + Sr ,, + u, i , ^ i + « , i , ^)- J k,k ( 29b ) 

Subtracting the mass average of (29b) from (29b) itself, we obtain 

m + = <u i ~ u i l£. + p ' lj k,k “ < ^" lj k,k > (29c) 

Multiplying (29c) by T" and mass averaging, we obtain: 

<T "uf" > + (T"^ - T " uj" )C j = <p- 1 T" J k k > - T TT </?- 1 J k;k > +.. (29d) 

where by ...(higher order terms) we mean all the terms that entail correlations higher than 
the second-order terms under consideration. For example, if we neglect the h o., we must 
also neglect uF in (29d): in fact, because of the second relation in Eqs.(5b), uF is already a 
second order variable As for the equation for T", we employ Eqs (27) and (32) of Canuto 
(1993; with obvious change in notation) to which we must add the last term in (15a). 
Keeping only the largest terms, we have 

gT"=-»;T r (uv T ''-^. + ,r; kk 

25 R .. IT /SC" j. SC 3T"s , 29 .. 

_ ~S Rc p *c U i ( 3* k W k W k > (29e> 

Once we multiply by c" and mass average, we obtain: 
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<c"DT">=-5X T ,i + X< c ''C2T" > + h o 


(290 


Adding Eq. (29d) to (290, we obtain 


g t TV = - *,T , - c p l HjC j + X<c"|>’> + <t> ‘T"J k k > - T-<p-‘J k>k > 


(29S) 


The last term becomes 


*cK T 


(29h) 


whereas 


A< c "^ 2 t " > + <f‘T"J kJ[ > = X<c"£2T"> + ^ c <T”£sc"> (29.) 

has a form analogous to (24c) and will be treated in similar fashion giving rise to the two 
most important terms 




(29j) 


Thus, finally: 

E.W t = -$T.-c; 1 H C. 
.Dt 1 p 1 ? ] 


r-JT'c" + Kx+X r )^ 0 T^- XrWT^-o^") 


dx? 


Q> T' 


PC 

<2x? 

(29k) 


VII. Non-Local Model 

To simplify the use of the equations we have derived, we list them below, beginning 


with the equations for the mean quantities: 


Large scale flow, U-: 

*$t u r - !xj( pi ij + * ,R ij ) 

(30a) 

Mean temperature T: 



- DT 

pc pTn ~ 

- (?Hi + rf)j + + “(Vp H i - ° c $ i) p ,i 

(30b) 
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where u is given by (15b) 

Mean concentration C: 

Reynolds Stresses /ra“uj = pR^: 




4 tVV = VVviN 


1J 


where 


- A ij ■ * + R jk u i,k) 

B ,j = 'ISVi - “c $ il p ,j + <H> 

The pressure-velocity correlation tt- is discussed in Appendix A 

Turbulent kinetic energy K=±R~: 

D 


(30c) 

(31a) 

(31b) 

(31c) 


(32) 


fcK+tyK) = - RjjUij - r l^Cp'H, - afpj - « 

In both (31a) and (32) we have not included the dilation term p 1 ^ 

Convective flux, CpPu"T" = pH. = CppJj: 

Dt H i + c p D f (H ,> = - Wr H j u u - y-’iH* - 

“ r pif H i + ?( i/ +^) H i ) kk 

Temperature fluctuations, ±pT" 2 = p\k: 

CT + ?' l Df(*) = -c p , H 1 T 1 -2 r ^ + ^ kk 

Concentration variance, ^pc" 2 = p$ 

g t *+D { («) = -* i C i -2r->* 

Concentration flux, pc"u"j = p$-: 

m*. +D f<*i) = - R ij C ,j - *j U i,j - ?-K^-2a c *)P,j - ^pj*i + *JC c <l+*'/X c )*i >kk 

(36) 

Temperature-concentration correlation, T"c": 

g t T^ + D { = - # iT - c p ‘H,C . - ; c p + - 


(33) 

(34) 

(35) 


- x c (o T T m -a r FT”)^ 




(37) 
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The time scales r pc , r c , r. 


Tq will be discussed below. 


VIII. Diffusivities. The K — e Model 

A widely used turbulence models is the non-local K— c model in which both K and t 
are treated non-locally while all the remaining turbulence variables are treated locally. The 
equations for the mean variables are unchanged We have two non-local equations: 

Kinetic energy K: 

Bt K + D f = -Vij +g Wi- Q c*i>-< < 38 > 

Dissipation rate e: 

Bt £ + D f = - c sVi,,+ C 1*¥ Vi-Vi^'-V 214 ' 1 (39) 

while the other turbulence variables are given by the local expressions 
Convective flux, F? = c^pu^'T" = c^J-: 

T P i> J i = - V J- ~ < 4 °) 

Temperature fluctuations, ±pT ' 12 = 

* = -iV, T ,i < 41 > 

Concentration variance, ±pc " 2 = 

$ = -ir c $jC. (42) 

Concentration flux, pc"u"j = 


- - R .j c ,j -*jU|J : 2» c *) p ,i 


(43) 


Temperature-concentration correlation, T"c'': 


r ., T n F= _ $i T.- J .c, (44) 

Reynolds stresses (Appendix A): 

wH (45) 

2r pv lb ij = -T5 KS ij - o-pp^ij - d-p 2 ) z ij + < 46 > 

Solving Eqs.(42)-(44), we obtain: 


dC 



where- 


r P ; d ik = R ik + < 47t>) 
"pc’ij = U i,j - 8 A i(-“ T T c9 T ,j + t c°<M) (47c) 

Aj = - (g?)-'§ (47d) 


Equation (47a) begins to acquire a familiar form but to obtain an explicit form for we 
must apply the Hamilton-Cayley theorem. The result is: 


•e* 

M • 

II 

1 

a 

^XjO 

(48a) 

where the turbulent diffusivity tensor D-j is given by: 


D ij = A (\ { ,k + Vik + ’imWA, 

(48b) 

with 


A = 1+L -L , A = - 1-L , A = (A + L Y 1 

0 12 1 1 0 3' 

(48c) 

L rV 2L 2 =- L l + Vji- 

6L 3= L ? + 2 ^im^mk% ” 3L , Vji 

(48d) 

From Eqs (42) and (45), we then obtain the expressions for the concentration variance $ 

and the T"c" correlation: 


fl._j._rv dCdC 
2 r c D ij 35E35T 

TVT =- r cA- T ,k D k i )i i 

(49a) 

(49b) 

Analogously, using Eqs (41), and (49b) into Eq (40), we obtain 
convective flux Jj which is structurally similar to (47a-c), namely, 

an expression for the 

^ik + ^ik) J k = " c ik T ,k 

(50a) 

where. 


r p0 c ik = R ik + °cS r c^i D kj$£. 

(50b) 

r p^ij = U i,j “S A i(~ r ^ a T T ,j + T c6 a M) 

(50c) 

Using the Hamilton-Cayley theorem, we can solve (50a). The convective flux is given by. 

J i “ " ^ik T ,k 

(51a) 
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The turbulent conductivity tensor ^ has the following form: 


X|j = B <Vik + Vik + Vmk> c kj 


with 


(51b) 


B = l+M -M , B = - 1-M , B = (B + M ) 

0 12 1 1 0 3 

M r^ 2M=-MU Vji , 

6M 3 = M’ + 2(‘ im % k ^ 1 , i - 3M 1 ( J IJ /‘j 1 


-1 


(51c) 

(51d) 


Finally, Eqs. (48a) for 4> ; and (51a) for Jj must be substituted in Eqs.(A.6-7) so as to 
obtain the tensor Bjj, Eq (A.5) Once that is done, the result is substituted in (46) and the 
Reynolds stresses R,. j can then be obtained in terms of the gradients of the mean variables 
The solution of (46) entails a system of algebraic equations. We recall that there are only 
five independent components of R since the kinetic energy K satisfies a separate 
differential equation (32). 


IX. No Mean Shear 

The ID case is particularly interesting since it allows a completely analytical solution 
of the problem. Using Eqs (48a) and (51a) one obtains the heat and concentration (salt) 
fluxes as 

= = (52) 

/>4> 3 = pw"s" = -pK c |§ (53) 

The turbulent diffusivities are given by the expressions: 

’ Kh^A’ K s=A < 54 > 

where the turbulent viscosity is given by 

v = rw 2 (55a) 

Ajj = 7^(1+ rpi + tt^xR )D _1 (55b) 

A g = 7^(1+ //x — tt 7r x)D'* (55c) 

D = (l+7pc)(l+//x) + (55d) 

1= p ), »= (55e ) 
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where we have introduced the following dimensionless functions 


x = r^N 2 

X i,2, 3) 4,5 = ( V V r 


pfl’ t 9> 


r-l 


(56a) 

(56b) 


where N£ has been defined earlier, Eq.(2b). Eqs.(54) are still not the final form since they 
depend on two unknown variables and x which we must express in terms of the large 
scale variables. To compute v T , we need an expression for w 2 For that, we use the equation 
for the Reynolds stresses, (46), (A.5)--(A.9). We obtain: 

"t = M 1 + 15< v A sV xl '‘ ( 56c > 

Next, we need an equation for x. We shall take the local limit of the kinetic energy 
equation (38) which reads 

^gVa-^sV^htAsVV (57a) 

Substituting (57a) into (56c), we obtain the equation for x 


15 


which changes (56d) to: 


*< A sV A h> - 7 


7 2 28 K' 

= T5 ej2 = T3T 


(57b) 


(57c) 


(57d) 


Thus, is expressed in terms of K and t. Finally, using the expressions for A^ c , Eq (57c) 
becomes 

A(x)x 2 + B(x)x - y 5 = 0 
where A(x) and B(x), which can depend on x (see below), are given by 

A=r( f i-i j r_ i )R () - x + 1 *^) (57e) 

B = (57f) 

Thus, x is expressed entirely as a function of R^. Finally, we have 

K _ 28 K 2 a K _ 28 K 2 

K h“ 15 7 A h’ K s _ 15 7 A s (58a ^ 

where we still have to determine K and t which in principle are solutions of the two 

dynamic equation (38) and (39). Eq.(38) has already been used in the local form, that is, 

Eq (57a) Eq.(39) for e has not yet been used and it can be taken to be local or not Below, 
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we give the solution corresponding to the case where (39) is taken to be local which means 

e = A- J K 3 / 2 (58b) 

where A is a mixing length; the specification of A is the price that one has to pay for not 
solving Eq.(39). From the definition of x, Eq.(56a), and the definition r=2Ke' 1 , we obtain, 
using (58b), 


K = 4A 2 N£x' 1 


and thus the final expressions for the diffusivities c follow from Eqs.(58a): 


K h= 13 A2 ^ A h- K s= 13 A2 (?)* A s 


N 2 


(58c) 


(58d) 


Thus, the problem is completely solved analytically. In fact, both diffusivities are now 
expressed in terms of the gradients ffT/dz and dS/dz Clearly, when N^<0 (corresponding 
to unstable stratification), x must be taken as the negative solution of (57e) since K is 
always positive, Eq (58c) 

In addition to the heat and salt turbulent diffusivities, it is also useful to introduce a 
mass diffusivity K^. We begin by using Eqs.(2b) and (22e) to rewrite (38) as 

§T + D f( K ) = K m N u - < 59a > 


where 


= p W (59b) 

is the "mass flux". Quantifying the strength of shear by the dimensionless parameter T 
(Hamilton et ah, 1989) 


r = knVM 

m u 

(59c) 

we have in the stationary and local case 


p' w" = pg A Te 

(59d) 

If we further write the mass flux as 


p'w" = - K^| 

(59e) 

the " turbulent mass diffusivity " becomes (Schmitt, 1994, Eq.7): 


K, = r ^ N’ = N£(i-R p ) 

(59f) 
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On the other hand, use of Eqs.(22e), (52) and (53) gives 


S^=Ng(K h -K c R p ) 


(59g) 


Using Eqs.(59e) and (2b), we obtain the following expression for K p in terms of and K c - 

K „ = ( R h - K cV (1 -V' (60) 

In the absence of mean shear, r= — 1, Eq (59d) shows that the mass flux is downward 

pV'cO (61) 


Thus we have, using (2d-g): 
SF, Stable, R p < 1 : 

Eqs.(59f) and (60) imply that: 


SF, Unstable R p >l 

Eqs (59f) and (60) imply that - 


v°- K i c <R p <1 


(62a) 


K p >0 ’ R p > KJ R p >l 


(62b) 


DC, Stable R p > 1 : 

Eqs (59f) and (60) imply that: 


K p <0, 


< > V 1 


(62c) 


The requirement of dynamical stability N 2 >0 sets the lower limit for R p while the 
requirement of turbulent mixing sets the upper limit of R . This is a natural result since 
transgressing the upper limit would mean that dS/di, which acts like sink, is too strong for 
turbulent mixing to survive. 

DC, Unstable, R p < 1 : 

Eqs.(59f) and (60) imply that: 


K p >0, 


> R . R <1 
P P 




(62d) 


X. Qualitative Results. No Shear 


23 



Before presenting the numerical solutions of the model, we present some qualitative 
results. Using the definitions of g , Eqs.(58d) and (57c), we derive the relations. 

w: = R p-^A; (63a) 

s r s 

In dtfjusive-convection , x<0 and since in the stable case R^>1, we conclude that 

K h >K g (63b) 

in accordance with the measurements (Kelley, 1984). In salt fingers , we write (57c) as. 

XT =R /5 1 ( 1+ 7 5 ^J (64a) 


Since x>0 and R^<1, it follows that 


K 0 > Kv 
s h 


(64b) 


in agreement with the measurements (Hamilton et al , 1989, fig. 2). Furthermore, in 

diffusive -convection, the flux ratio 

o$ K c ,, 

R = _£3 = s R R ( R 15 x -i A -i)-i (65a) 

F a T J s K h p p K p 7 s ' v ' 


is predicted to be (x<0) 


R p <l (65b) 

m agreement with the data (Kelley 1990, fig. 2). Similarly, in salt fingers we derive that the 
flux ratio: 

r f = aV = rV = (i+^x'^h 1 )* 1 ( 66a ) 

c 3 s K 

is predicted to be (x>0) 

R F <1 (66b) 

in accord with the measurements (Turner, 1967, fig.4; Schmitt, 1979, fig. 4; McDougall and 
Taylor, 1984, fig.4; Taylor and Buchens, 1989, fig.6; Ozgokmen et al. 1998, fig.13). 


XI. Salt— Fingers and Diffusive-Convection. The effect of Shear . 

Here, we present the analytic solutions for the turbulent diffusivities of momentum, 
heat and salt in the presence of the three gradients VU VT, and VC which we take in the 
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form given of Eqs.(73a,b). K and e can be treated either locally or not. 
introduce the following dimensionless variables: 

It is convenient to 

n , = - v 3 ga T 72 $£ i 

(67a) 

c i e i ^“slj 

(67b) 

% = ^ K 'V T pv J i 

(67c) 

</>=/? K’*a g r $. 

*1 m 5 c® pv 1 

(6’d) 


(67e) 

eqs.(46), (47a-d) and (50a-c) become: 


Reynolds stresses: 


v K 'V^j 

(68a) 

^ij = - T5 £ ij “ - (l-P^y + V T lj 

(68b) 

*ij e A.^ + A,Pj - jA k ^j» k 

(68c) 

T ij = A i^j + Vi-iWk 

(68d) 

Concentration Flux: 


!< ik + ’ik^k = “ l p 4 (a ik + Plk'l + P s A i a k | l C k 

(69a) 

"ij = P 3 C i.j “ A iV”j + C j> 

(69b) 

Temperature Flux: 


(<i ik + "ik^k = p 6 < a ik + 3 { ik “ p 7 Vk )n k 

(70a) 

"ij = p 8 °i.j - A i( p 9 p j + p 10 C j) 

(70b) 


where, 
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( 72 ) 


P = 7T 7T 7 1" 2 , 
10 2 4 3 


D = 7T 7 T' 1 
11 1 3 


If we take. 


and thus- 


1 (T,C) - «, 3 |(T,C), u = [U(z), V(2), 0] 


(73a) 


V* 


rO 0 dt}f dz-i 

0 0 dV/dz 

dtS/dz dV/dz 0 


V..=i 

’ v y 2 


0 0 dX}/dz - | 

0 0 / dz 

— dt)/dz- dV/dz 0 


(73b) 


we can give a complete algebraic solution of the algebraic set of equations (68)— (70). Since 
we are dealing with only one component of the vectors n^, Cj, we simplify the notation to. 

(74a) 
(74b) 


id T 


n =n= - n go T r 2 ^, n = 7r 7 t 

3 o T <7Z’ 0 2 3 

c =c=c gT 2 a ^, c = n 2 
3 o to cdz’ o 3 

The dimensionless shear is given by: 

V N u> 2 ’ K = 

If we introduce the simplifying notation 

7=7-1 /) — r P I * 


w=/rVu z , 0=T' 


we obtain the following results: 

Momentum Flux. 


(74c) 

(74d) 



uw - ■ K m = 2 7 S m 

(75a) 

Heat flux. 

^ = " K hl?’ K h = 2 7 S h 

(75b) 

Salt Flux: 

^ ,= - K sI’ K s= 2 T S s 

(75c) 

The dimensionless structure functions S , , see Eqs.( 3 d), are given by 



S S, = U A.D- 1 , S = Yc7r A D * 1 
m 75 m ’ n 15 4 n 5 s 15 1 s 

(76a) 


=12 + an 2 + anc + ac 2 + an + ac 

m 1 2 3 4 5 

(76b) 


= (1 + b^c + b 2 n)(60 4- b g y + b^c + b^n) 

(76c) 


A c = (1 + b^c + b^n)(60 + b g y + b^c + b^n) 

(76d) 
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D = 24 + d yn 2 + d ync + d yc 2 + d n 3 + d n 2 c + d nc 2 + d c 3 + 

1 2 3 4 5 6 7 

+ d yn + d yc + d n 2 -f d nc + d c 2 + d y 4- d n + d c (76e) 

As one can see, the dimensionless functions A's and D depend on the gradients of the mean 
temperature, concentration and mean velocity represented by n, c and y. The functions a^, 
and d^ (Appendix C) depend on the time scales r c , r pc , etc which in turn depend on the 
Peclet numbers. For large Peclet numbers, a , b and d become constant, Appendix C. As 
before, the variables K and e are in principle solutions of Eqs.(38) and (39). The 
(superficial) algebraic complexity of the functions A's is a small price to pay when one 
considers that the above equations are the solution of a fully turbulent problem in the 
presence of three external fields, T, U and C. It is indeed quite surprising that such a 
complex problem could be expressed via a set of algebraic relations. 

In the case of a local model, Eq.(38) becomes: 


-R i j U i J + 6a T Vi-S a s $ i = t 

Using the definition of the K's and that of ip given in Eq.(74c), we have 


(77a) 

(77b) 


Once we substitute the functions S m ^ , Eq (77c) yields the function: 

ip=iP{ Ri,R p ) (78a) 

We recall that in the functions A^ u „ we must substitute: 

m,h,s 

n = -| 5 n^Ri(l-Rp-> (78b) 

c = l-Rp- 1 (78c) 

We shall exhibit the turbulent diffusivities K m , and K g in units of A 2 N u (for 

different values of R ) vs.Ri which we recall is defined as follows: 

P 

Ri = ( S a T !£) Nu ' 2 ( 78d ) 

which helps us differentiate between stable and unstable stratification. 


XII. The RNG method to determine the time scales r pc , t q , Tq 

To make the above equations predictive, one must know the dissipation time scales of 
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the different turbulent variables, namely r , t ^ r c , r Tq. Not surprisingly, this is one 
of the most difficult problems since one-point closure models, like the one we have used, 
are unable to provide them. In most engineering and geophysical applications (e.g., the MY 
model), it was always assumed that 

\ = ( v r c 0 ’ T e V T e > r '~ consUnt < 79a > 

However, on physical grounds, it is only possible to say that 

T P r V T r T c (79b) 

while remains to be determined. Since in principle, one may want to consider regimes 
in which the Peclet number of both the temperature and salinity fields are not excessively 
lager than unity (Pe~l correspond to low levels of turbulence), we adopt expressions for 
that were previously depermined: theye are 

r pc )r‘ = aPe(l+bPe)-‘ 

(79c) 


47 r 2 a =1, 5a(l+cr ^ 1 )' 1 

( r 0 , r c ) 7 " 1 = aPe ( 1 + aPecr t 1 )' 1 

7 A =4 


(79d) 


W‘ = &Pe «( 1 + b ‘'ti 

7 T^aE 4(1+Pe^/Pe c ) _1 

4b = 15a(l+a t ^/(7 tc ) (79e) 

We have used only one symbol for both Pe and a ^ but clearly in each specific case one must 
insert the corresponding Pe and <r t , where: 

Pe 0 ,c = m 7 ^c 1 ) (79f) 

where Xq z are the molecular diffusivities of the two fields. The turbulent Prandtl numbers 

<r t are functions of the corresponding Pe and the RNG method gives the following result 

(Canuto and Dubovikov, 1996) 

vi' = 1 + !"V e ' 1(|1 + l? Pe(<r i 1+7 ; 1) ' r (79g) 

The constants 7 are given by: 

1?2 
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2? t = (T 2 + 47 )*— 7 , 7 2 =7j+7i r= 7 1 /7 2 , 7=0.3 (79h) 

XIII. Numerical Results 

In Figs.1-3 we plot K m ^ g vs. Ri (defined in Eq.78d) for different (defined in 
Eq.2a). The panels are characterized by the symbols SF (salt fingers), DC (diffusive 
convection), DS (doubly stable) and DU (doubly unstable) defined in Eqs (2d-g). Consider 
the case of salt-fingers in Fig.la. At a fixed Ri, the diffusivity increases as R^ increases 
which is physically understandable since the instability is generated by salt and thus the 
larger the source, the larger the diffusivity. Next, consider the dependence on Ri. We notice 
that the smaller the shear Ri-*oo, the larger are the K's, which at first may seem 
paradoxical: since both salt and shear contribute to the instability, their effect should add 
up What we find is that the larger the shear, the smaller the diffusivity, which implies 
that shear and salt-fingers work in opposite directions. It is in fact known (Linden 1971, 
1974b, Kunze, 1990) that shear has the tendency to disrupt the fingers transport process In 
the case of DS and DU, Fig.lb, R^ is negative, see Eqs.(2f-g). Quite understandably, the 
former case (right panel) corresponds to the lowest diffusivity because of the large stability 
introduced by both salt and temperature. The only source of instability is shear and thus 
turbulent mixing dies when stratification is too strong. In the DU sade, the opposite occurs 
in the sense that both T and S are unstable and the resulting diffusivities are the largest. 
The same considerations hold true for and K g which are shown in Figs. 2, 3. Consider 
now the DC case, Fig.la. At a given Ri, the diffusivity decreases as R^ increases, the 
opposite of the SF case. This is in accordance with the fact that in this case salt acts as a 
sink of turbulent mixing (which is caused by an unstable temperature gradient), and thus, 
the stronger the sink, the lower the level of turbulence, a circumstance that is reflected in 
the decrease of the diffusivity. As for the effect of shear, we notice that here too, the 
smaller the shear (large Ri), the larger the diffusivities which implies that shear prevents 
the mixing caused by the temperature instability. However, this is not true in general, the 
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curves first decrease with increasing Ri, which indicates that for moderate Ri shear helps 
mixing, as one would expect, but the trend does not continue since the curves change 
curvature. However, there is saturation phenomenon which does not occur in the SF case 
At large R^ (large sink), the help in mixing from shear saturates. Finally, the lowest three 
curves correspond to a stable situation, while the second and third correspond to an 
unstable situation. In Figs. 2b and 3b we present and K g in the DU (doubly unstable) 
and DS (doubly stable) cases. Quite naturally, in the altter case the diffusivities are the 
lowest. In Figs.4-6 we plot the ratios K m /K g and K^/Kg which show quite 

clearly that the diffusivities are indeed different among themselves. In Figs. 7-10 we exhibit 
the turbulent mass diffusivity K defined in Eq (59e) and given in terms of K. by 

P 

Eq (60). In Fig 11 we plot T defined in Eq.(59f). Schmitt (1994) "measured" values of 
T=0 18-0.25 are indeed predicted by the model for the case of salt fingers (upper right 
panel) for quite a range of Ri but the precise value depends on R^. 

The length scale A is determined using the Deardorff-Blackadar formula: 

A = 2‘ 3 / 2 B l, B =24.7, 

1 1 

<=nrtn(*ft,<) (80a) 

( = at (t +ra) J , t = 0.17H (80b) 

where iq 2 =K is the turbulent kinetic energy, N is the Brunt-Vaisala frequency, k= 0 4 is 
the von Karman constant and H is the mixed layer depth. When used within the NCAR 
CSM Ocean Model, H is determined as the depth where the buoyancy difference 

g[p(H)- p(surface)]p(H)‘ 1 = 3 lO'W 2 (80c) 


XIV. Ocean GCM 

We tested the new vertical diffusivities in a global ocean general circulation model, 
the NCAR CSM Ocean Model produced by the University Corporation for Research, 
National Center for Atmospheric Research, Climate and Global Dynamics Division. They 
developed their model by modifying the MOM 1.1 GFDL code (NCAR CSM Ocean Model 
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Technical Note, The NCAR CSM Ocean Model, by the NCAR Oceanography section). We 
employed the stand-alone 3°x3° configuration of the model detailed in their technical note 
with the default parameter values. It has 3.6° spacing in longitude and a variable spacing in 
latitude increasing from 1.8° at the equator to 3.4° at 17° N, S and then decreasing back to 
1.8° for 60° N, S and poleward. There are 25 levels of increasing thickness in the vertical, 
with the surface level 6 meters thick The option for the GM mesoscale eddy 
parameterization was enabled. Bulk forcing with a seasonal cycle plus a 1/2 year timescale 
restoring condition on the salinity is used, except under sea-ice where there is strong 
restoring. This configuration corresponds to the case B-K described in Large et al (1997). 
It should be noted, however, that for determination of the length scale in our turbulence 
model we used the program's definition for mixed layer depth (a buoyancy difference from 
the surface of 3 10‘ 4 ms' 2 ), which is different from that graphed as a diagnostic in Fig. 5 of 
Large et al.(1997) We initialized our runs with annually averaged Levitus data and ran for 
126 momentum years. As in Large et al. (1997) a 3504sec timestep for momentum is used, 
while for the first 96 momentum years the tracers are accelerated by a factor increasing 
from 10 at the surface to 100 for the deep ocean. We then set all timesteps equal for the 
remaining 30 years as they did. 

First, we ran the NCAR program as is, with the option for the KPP mixing enabled, 
producing the KPP data presented in the figures below. Then, in place of the KPP module, 
we inserted a module which uses our new model for the diffusivities for momentum and 
heat with the salt diffusivity set equal to that of heat. To save computing time, we 
constructed tables of the dimensionless functions S m ^ and of the dimensionless variable y 
(obtained from solving Eq.68e), 

y = i nr = iA*) 2 (81a) 

vs. Ri. Then, for each point in space and time these were interpolated to the local Ri. The 
diffusivities K m ^ were written in terms of (81a) as 
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K m,h/ ft - * B r S m ,h 


(81b) 


XV. Below the Mixed Layer 

Below the ocean mixed layer, the external wind-generated shear is too small to 
generate turbulent mixing and yet, even in regions where both the temperature and the 
salinity gradients are stably stratified, it is usual to assume "background diffusivities" for 
viscosity, heat and salt diffusivity which are believed to be caused by internal wave 
breaking (Large et al., 1997). In our case, when we assumed Kj 1 =K g , we followed the same 
practice. It would be preferable not to do so but rather model the physical processes 
causing this background mixing. Our main assumption is that the turbulence model has 
given us the correct functional dependence of the K m ^ g on Ri and and that such 
diffusivities can thus be used below the ML. Since all the arguments discussed below, are 
valid for any of the three K's, we shall use only the generic symbol K and write succinctly 

K = K(Ri,R ) (82a) 

The key problem is how to define and thus compute Ri. Here, we shall make us of the 
measured data (Gargett et al., 1981) on the vertical shear generated by the wave breaking 
phenomenon. By integrating over all wavenumbers one can compute the shear due to 
internal waves, S^. One can then form a corresponding Ri^ as follows: 

Ri wb= N2 /S£b ( 82b ) 

where 

n j = ( 82 °) 

Gargett et. al. (1981, sec. 5) confirmed earlier arguments by Munk (1966) that Ri wl3 ~l To 
those argument, we would like to add the following consideration. As the value of Ri cr , 
above which there is no longer turbulent mixing, computed from our model is 0(1), if Ri w ^ 
were >>1, there would be no turbulence generated by the internal waves at all. On the 
other hand if Ri w ^ were <<1, there would be a very strong turbulence producing a 
viscosity sufficient to damp out the waves themselves. The wave-generated turbulence is 
thus self-limiting Since the turbulence model gives a precise value for Ri cr , while the 
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above argument only tells us that Ri w k~0(l), we shall write: 


Ri , =cRi 
wb cr 


(82d) 


where c is a constant reasonably close to unity. We have found that c=0.88 gives a 
diffusivity close that measured by Ledwell et al. (1993). Since in the local model, the K's 
are also proportional to the length scale A or 1 . Below the mixed layer, we thus need an 
analogous £ wb' We shall use the same formal expressions (80a, b) but with different ^(wb) 
which we compute as follows. Assuming a Kolmogorov spectrum at wavenumbers upward 


of a breakpoint and integrating, we obtain: 


»b) = (3 Ko) 3 / ! (B k 4 )-' 


(82e) 


where Ko=1.6 is the Kolmogorov constant. We identify k Q with the best value of Gargett 
et al. (1981) for the break in slope of the observed spectrum of internal waves, namely 

k^ = j-Q 2 tt radians/meter (82f) 

Thus, f Q (wb) is known and so is W Similarly, y w ^ is obtained by solving the 
production=dissipation, Eq (77b). Thus, the complete wave-breaking expressions for the 
three diffusivities are: 

K m,h,s( wb > = * B /( wb >Wwl, S m,h,s( Ri wb' V < 82 8) 

We add together the diffusivities calculated using the shear resolved in the ocean model 

and the background diffusivities, ensuring continuity in the transition between regions 

where external excited shear dominates and those where the internal wave shear does We 

thus take the total diffusivities to be: 

K m,h.s = K m,h,s( Ri 'V + K m,h,s< wb ' V < 82h > 

In the statically unstable case (Ri< 0), we set K , (wb)=0. The very large mixing due 

to convective instability makes the background irrelevant in this situation in any case. 


XVI. O-GCM results 

Using the model for the K's extended all the way to the bottom of the ocean, we 
obtain the results presented in Figs.12-23. In each case, we compare the results with 
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Levitus (1994) data, with the KPP model (K^=K g ) for which we have rerun the code and 

with our model with K=Kv. In Figs. 24-32 we plot K , (cmV 1 ) vs. depth (meters) 

s n 

at different locations. As expected, the K's are small below the mixed layer where they can 

reach very high values, as we explicitly show in Figs. 30-32. In case of the Canary Islands, 

Fig.29, the diffusivity of a truly passive scalar (and thus strictly not K , ) was 

rn,n,s,p 

measured by Ledwell et al. (1993) to have a value of 0.11±0.02 cmV 1 . Finally, in Fig.33 we 
present the polar heat transport. As already discussed in the work of MHG, global 
properties are not strongly affected by double diffusion phenomena. 

XVII. Conclusions. 

Considering the importance of double diffusion phenomena in oceanography (Schmitt, 
1994, Zhang et al., 1998; Merrifield et al., 1999), we believe we have made a quantitative 
step by presenting a new formalism. The resiliency of the new approach is demonstrated by 
the fact that it can encompass salt-fingers, diffusive-convection, doubly stable and doubly 
unstable gradients The whole formalism was developed so as to include shear which, 
though of different origin at differernt depths, is always present. Within the mixed layer, it 
is mainly due to external wind gradients while in the ocean interior is believed to be mainly 
due to wave breaking phenomena. 

Clearly, the model with salinity would have lost much of its attractiveness if we could 
use it only in the ML and if we had to parameterize the physical processes below the ML 
with adjustable background diffusivities as done thus far. We have suggested that below 
the ML the functional dependence of the three diffusivities on the two stability parameters 
and Ri is still the one given by the turbulence model since the latter does not depend on 
any specific form of the shear entering the Richardson number Ri. As dicvussed in XIV, we 
have used the data on vertical shear measured by Gargett et al. (1981). 

The final model comes in more than one flavor depending on whether on uses local or 
non-local models. Logically, the first model we have treid is the local one since the whole 
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problem can be solved analytically. The expressions for the turbulent diffusivities are 
algebraic. The whole turbulence problem is reduced to the solution of a cubic equation 

The problem is however far from solved. Of particular significance is the role played 
by the salinity-temperature correlation. If we were to assume that 


TV = (T^T 2 )* 


(83a) 


and that 


r c r * T c = * T t 


(83b) 


as one may be tempted to do, one would obtain that the two fields are indistinguishable 
and this implies that 


K g =K h (83c) 

contrary to what is observed. Fortunately, the present model dos not require either of 
Eqs (83a, b) but once they are imposed, (83c) follows. These and similar questions will be 
the subject of future studies. 
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Figure caption 


Fig. la Momentum diffusivity K m in units of A 2 Nu, see Eqs.(74c) and (58b) vs. Ri defined 
in Eq.(78d) for different values of the Turner number defined in Eq.(2a). The label DC 
and SF are defined in Eqs.(2d)-(2e). 

Fig. lb Same as in Fig.la for the DU and DS cases, Eqs.(2d)-(2e). 

Fig. 2a. Heat diffusivity vs. Ri for different R^. Salt-fingers and diffusive-convection 
Fig.2b Same as in Fig.2a for the DU and DS cases. 

Fig. 3a Salt diffusivity K vs. Ri for the SF and DC cases 

o 

Fig. 3b Same as in Fig.3b for the DU and DS cases 

Fig 4 The turbulent Prandtl number K m /Kj i vs. Ri. The heavy line corresponds to the 
laboratory data discussed in paper II, Figs. 3,4. 

Fig 5. The ratio of momentum to salt diffusivity vs. Ri for different R^ 

Fig. 6 The ratio of heat to salt diffusivity vs. Ri for the DC and SF cases for different R^ 

Fig 7 The mass flux diffusivity defined in Eqs.(59e) and (60) vs. Ri for different R^ for 
the DC and SF cases 

Fig. 8 The ratio K^/K m vs. Ri for different R^ for the DC and SF cases. 

Fig. 9 Same as in Fig.8 for the ratio K^/K^ 

Fig 10 Same as in Fig.8 for the ratio K /K 

p s 

Fig 11 The efficiency parameter T defined in Eqs.(59c) and (59f) vs. Ri for different values 
of R . A value r=0.18-0.25 (Schmitt, 1994) is indeed predicted by the model for the 
salt-finger case. 

Fig. 12 The resulting global ocean temperature using the O-GCM discussed in XIV with 
the background diffusivities computed following the new procedure discussed in XIV above. 
The Levitus (1994) data are the solid line. We have also run the 0— GCM code with the 
KPP model (K^K^) and the results are indicated by diamonds. The results with our new 
model with Kg=K^ are shown by squares while the full model with K g #K^ are indicated by 
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asterisks. 

Fig. 13 Same as Fig.12 for the global salinity 

Fig. 14. Same as Fig.12 for the Artie ocean 

Fig.15. Same as Fig 13 for the Artie ocean 

Fig. 16 Same as Fig.12 for the Atlantic ocean 

Fig. 17. Same as Fig. 13 for the Atlantic ocean 

Fig. 18. Same as Fig.12 for the Pacific ocean 

Fig. 19. Same as Fig. 13 for Pacific ocean 

Fig. 20. Same as Fig.12 for the Indian ocean 

Fig. 21. Same as Fig. 13 for the Indian ocean 

Fig. 22. Same as Fig.12 for the Southern ocean 

Fig 23 Same as Fig. 13 for the Southern ocean 

Fig.24. The foir diffusivities K , (cmV 1 ) for the Papa staion 

Fig. 25 Same as in Fig.24 for the Artie ocean 

Fig. 26 Same as in Fig.24 for the Canary Islands. 

Fig 27 Same as in Fig 24 but for the first 1km 
Fig 28 Same as in Fig. 25 for the first 1km 

Fig. 29 Same as in Fig 26 for the first 1km. Ledwell et al. (1993) value of 0.11±0.02 cm 2 s 4 
(see, however, the discussion in the main text) 

Fig. 30. Same as in Fig.27 for the first 40m 
Fig.31 Same as in Fig.28 for the first 40m 
Fig.32 Same as in Fig.29 for the first 60m. 

Fig. 33 Polar heat transport vs. latitude for three different models. 
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Appendix A: Reynolds stress equations 


Rather than employing Eq.(31a), we introduce the traceless tensor 


WsVkk =R ij-3 { ij K 

(Al) 

where K satisfies Eq. (38). We thus have: 


Dt b ij + D f(b) = - f K£ ij - fly- Zjj + By - 

(A.2) 

where the (traceless) tensors ft and Z representing shear and vorticity 

are defined as: 

f! ij = b ik E jk + b jk E ik “ hWu 

(A 2) 

Z ij= b ik V jk + b jk V ik 

(A.3) 

where and V- are shear and vorticity: 


V *< u ij + V- V ii = * (u i.j-V 

(A. 4) 

The new tensor B- is given by 


B ,j = 6(“T L ij - “c M ij) 

(A 5) 

VVj + Vrs^jVk 

(A 6) 


(A 7) 

We recall that 


A i =-&% 

(A 8) 

Finally, we have to treat the pressure— velocity tensor. Following the procedure described in 


(Canuto 1994), we take 

= 2 v‘ b o - K - p .v p * z ij + (A - 9) 

where the numerical constants p and 3 are given in the text. The time scale r is 
discussed in Appendix B. Finally, Eq. (A. 2) becomes 

m b ij + D f(») = - 2 V% "is KS ij - z ij + ^ B ij 

(A.10) 


Appendix B 

The (r , r Tq) vs. r relation is (Canuto and Dubovikov 1998): 
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t = 2Kc\ 


r pv ” 5 r 


for the T-field we have 


^ 6 - t 1 + 


T~ ~ 7^ Pe 0 f 1 + 7^*8°'$ 


n-i 


For the C— field we have: 


(B.l) 

(B 2) 
(B.3) 

(B 4) 
(B.5) 


4M 


r^ C ~ l^ ?e c f 1 + 

t ~ = 7i 2Pe C l 1 + 7i 2Pe c <r tcl‘ 1 
For the T-C correlation, we have: 

F* = ^(I+Pypep-'U + ^Pe^J, (l+«r t< / V )(l + Pe/Pe c )-'] 

The Peclet numbers Pe^ are defined as: 

Pe - 471-2 K V 1 ) 

*>c " 125 c X C J 

The turbulent Prandtl numbers <7^, are themselves functions of the corresponding Pe's 

and satisfy the general equation. Calling < 7 ^=£, we have 

o E 7 £+1 p 

^ —9-ri __1 / . _ir/< .0 T«- 1 \-l 


(B 6) 


(B.7) 


7 2 E = 1 + ^Pe->(7 2 -o)[(l + I^Pe ) J -1] 


(B.8) 


4* 

with 27 i =(t 2 + 47 ) 2 - 7 , 7 = 7^+7 and 7 = 0 . 3 . The Prandtl number o=vj\ is usually O(10‘ 8 ) 
and thus negligible. 

The Peclet number Pe can safely be taken much larger than unity in which case both 
(B.4) and (B.5) become constant. When also Pe^>>l, we have 


a t =0.72 


and thus: 


or 


T J T=T J T= 3 (1+< ’t‘t 1 > V T=r c/ T = <V W T= I5°t 


V /r= V /r = 0 0837 


(B.9) 

(B.10) 
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Tl> ese val Ues . ^ T c /t = 0t 2 

UnimP,ytht ^<^Z° 96 

P 1 ==Q 832 n 

^0, 0n P f°*. 

P ^ P C 00323 

p ^o.4rg 9} 6 422 

p ^,r K 


We ^o have 

Thus: 


P = 

11 


^m^O-ies, 


p . 

10 


* V 

: 0 -ll63 


: 0 0155 


P 2 in' '0-455 


a r 10494, . n 
: °-0l63 » 2 " 0-9239 

b ’ 4SZ ' 10 - 4 205, 
i^'O-ioog ^ s ^-3656 

‘."OUtUB; ' b C' 0,163 

Vw.** / '°' 9 ® 9 

7'"0.?5 58 

fj — . n - 


•> d 
2 

^ o 1042 

d ^ 

4 

" 0.3494 

d * 
6 

" °-05?2 

d ^ . 
8 

' l 0938 

d ^ 
10 

6 -227l 


d 4 0Q{;n ’ V *227j 

P 4 -°950, w „ 

d = i jc^ J2 ' 0 ' 1034 
Quantities n ^ d J ‘^'^S 

’ Wp// ^ 

elIas (90b c ) * 




(B.12) 

(B 13) 


(B-14) 


(B.15) 


(B 16) 
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n = 0.0691, c = 0.5184 
o o 

Appendix C 

The functions a,b,c and d entering Eqs.(88a)-(88e) are given by: 
a = p [I2p +8p -30p p -5p (p +3p )], 

1 ll 1 9 6 6 8 6 V im *2m' J ’ 

a = - 5[(p p +p p -2p p p )(p +3p m )+ 

2 u *4 9 6 11 4 6 7 A l m 2^ 

+8(p 4 p+2p 6 p ir 2p 5 p 6 )+12(p ii p 9+ p )0 p ii -p 5 p 6 p 7 ) 

-30(p 3 p 4 p 9 +p 6 p 8 p ir p 5 p 6 p 8 -p 3 p 5 p 6 ) 

a = p [8p +12p - 30p p - 5p (p ^+3p _)] 

3 10 l 4 11 * 3*4 4 Vi im *2m /J 

a = - p (8 - 30p - 5p - 15p „) - 12(p +p ) 

4 6 V 8 im *2m' '*9 i v 

a = - p (8 - 30p - 5p — 15p ) - 12(p +p ) 

5 4 V 3 im * 2 m 7 ' io n 

b = p p -p , b =-p , b = 15p 2 +2p ^ -5p 2 -6p 

i 4*7 n 2 li 3 2 m *im *im 

b = - 30p , b = - 30p , b = - p . b =p p -p 

4 * 4 5 6 6 10 7 6*7 *9 


(Bn) 


2m 


(Cl) 


(C 2) 


d = p [p 2 (p +6p ) + 2(p — 3p ip p -p 2 (p +2p )] 

l *ir 2 m '*6 *V v *im ^ 2 m /K 6 K 8 *im v *e *V J 

d = (p 2 -p 2 )(2p p p -p p -p p ) + 

2 v *im ^2m /v *4*7*6 *6*11 *4*V 

+ 2 ^m- 3 P2m)(P 4 PeP?-P„Pr P ,o P u ) + 

+ 2 (P 1 m- 3 P 2 m)(P3 P 4 P ! , + P«P 8 P 1 rP 4 P 6 P,P 8 -P3 P 4P6P7 ) 

d 3 = p .o |p 2 m(P4 +6p „) + 2 ( p ,m- 3 P 2 m) P 3P4 " Py p , +2p „M 

d r- 4 P 6 P u ( 2 P 6 + 3 P 9 ) 

d = 4p p p 2 (4+3p ) - 4p p (3p +2p ) - 4p p (3p +3p +2p +2p ) 

5 4 7 6 7 4 9 11 6 6 11 9 10 4 6 * 

d = 4p 2 p p (4+3p ) - 4p p (2p +3p ) - 8p p (p +p ) - 12p p (p +p ) 

6 4 6 7 7 4 9 4 11 4 6 10 11 10 11 4 6 

d = - 4p p (2p +3p ) 

7 4 10 4 11 

d = p 2 (2p +2p +p ) - p 2 (6p +6p +p ) - 2p p (p - 3p _,) 

8 *im v *9 *11 *V y 2 in K *9 *11 *V * 6 * 8 V *im * 2 m 7 

d = p 2 (2p +2p +p ) - p 2 (6p +6p +p ) - 2p p (p - 3p _J 
9 K im v *n *io *V * 2 m v K n F io *V * 3 *V*im F 2 m y 

d = 8p 2 + 4p (3p +7p ) + 24p p 

10 6 6 V 9 ll* 9 11 

d = - 8p p p (4+3p ) + 4p (4p +7p +3p ) + 4p (3p +7p ) + 24p (p -f p ) 

11 4 6 7 r *V 6 9 \\' *6 10 11 11 V 9 10 * 


d = 4p (7p +6p ) + 4p (2p +3p ), d = 6p 2 - 2p 2 

12 io 4 lr 4 V 4 lr 13 2 m *im 

d = - 24p -24p -28p , d = - 24p -24p -28p 

14 *9 11 6 15 11 10 4 


(C.3) 
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